COMMENT ⊗   VALID 00006 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	 Load MACHAR routine
C00011 00003	 Random number generators
C00012 00004	 sqrt Testing
C00022 00005	 atan Test Routine
C00039 00006	 Testing
C00040 ENDMK
C⊗;
;;; Load MACHAR routine
(eval-when (compile load) (fasload machar))

(declare (special *ibeta* *it* *irnd* *ngrd* *machep* *epsneg* *maxexp*
		  *negep* *iexp* *minexp* *eps* *xmin* *xmax*))

(declare (flonum *xmax* *xmin* *epsneg*)
	 (fixnum *ibeta* *it* *irnd* *machep* *minexp* *iexp* 
		 *negep* *ngrd* *maxexp*))

(declare (*expr machar))
;;; Random number generators

(declare (special *iy* *j*) (fixnum *iy* *j*)
	 (flonum (ran fixnum))
	 (flonum (randl flonum)))

(setq *iy* 100001.)

(defun ran (k)
       (setq *j* k)
       (setq *iy* (* *iy* 125.))
       (setq *iy* (- *iy* (* (// *iy* 2796203.) 2796203.)))
       (//$ (float *iy*) 2796203.0))

;;; Get value of E in here
(defun randl (x)
       (exp (*$ x (ran 1))))
;;; sqrt Testing
(declare (special *results*))
(setq *results* ())

(defun sqrt-test ()
       (declare (flonum beta sqbeta albeta ait a b xn c x1 r6 r7 x y z w)
		(fixnum n i k1 k2 k3))
       (machar)
       (let* ((beta (float *ibeta*))
	      (sqbeta (sqrt beta))
	      (albeta (log beta))(ait (float *it*))
	      (a (//$ 1.0 sqbeta))(b 1.0)
	      (n 8000.)(xn (float n))
	      (c 0.0)(k1 0)(k3 0)(x1 0.0)(r6 0.0)(r7 0.0)
	      (x 0.0)(y 0.0)(z 0.0)(w 0.0)(k2 0))
	     (do ((j 1 (1+ j)))
		 ((= j 3))
		 (setq c (log (//$ b a)))
		 (setq k1 0 k3 0 x1 0.0 r6 0.0 r7 0.0)
		 (do ((i 1 (1+ i)))
		     ((> i n))
		     (setq x (*$ a (randl c)))
		     (setq y (*$ x x))
		     (setq z (sqrt y))
 		     (setq w (//$ (-$ z x) x))
		     (cond ((> w 0.0) (setq k1 (1+ k1))))
		     (cond ((< w 0.0) (setq k3 (1+ k3))))
		     (setq w (abs w))
		     (cond ((< r6 w)
			    (setq r6 w)
			    (setq x1 x)))
		     (setq r7 (+$ r7 (*$ w w))))
		 (setq k2 (- (- n k1) k3))
		 (setq r7 (sqrt (//$ r7 xn)))
		 (push '(test of sqrt(x * x) - x) *results*)
		 (push 
		  `(,n random arguments were tested in the interval
		       (,a ,b)) *results*)
		 (push 
		  `(sqrt(x) was larger ,k1 times) *results*)
		 (push `(it agreed ,k2 times) *results*)
		 (push `(it was smaller ,k3 times) *results*)
		 (push `(there are ,*it* base ,*ibeta* significant digits in
				a floating-point number) *results*)
		 (setq w -999.0)
		 (cond ((not (zerop r6))
			(setq w (//$ (log (abs r6)) albeta))))
		 (push `(the maximum relative error  of ,r6
			      = ,*ibeta* ↑ ,w occurred for x =
			      ,x1) *results*)
		 (setq w (max (+$ ait w) 0.0))
		 (push `(the estimated loss of base ,*ibeta*
			      significant digits is
			      ,w) *results*)
		 (setq w -999.0)
		 (cond ((not (zerop r7))
			(setq w (//$ (log (abs r7)) albeta))))
		 (push `(the root mean square relative
			      error was ,r7 = ,*ibeta* ↑ ,w) *results*)
		 (setq w (max (+$ ait w) 0.0))
		 (push `(the estimated loss of base ,*ibeta*
			      significant digits is
			      ,w) *results*)
		 (setq a 1.0)
		 (setq b sqbeta))
	     (push `(test of special arguments) *results*)
	     (setq x *xmin*)
	     (setq y (sqrt x))
	     (push 
	      `(sqrt(*xmin*) = sqrt(,x) = ,y) *results*)
	     (setq x (-$ 1.0 *epsneg*))
	     (setq y (sqrt x))
	     (push `(sqrt (1.0 - *epsneg*) = sqrt(1.0 - ,*epsneg*) = ,y) *results*)
	     (setq x 1.0)
	     (setq y (sqrt x))
	     (push `(sqrt(1.0) = ,y) *results*)
	     (setq x (+$ 1.0 *eps*))
	     (setq y (sqrt x))
	     (push `(sqrt (1.0 + *eps*) = sqrt (1.0 + ,*eps*) = ,y) *results*)
	     (setq x *xmax*)
	     (setq y (sqrt x))
	     (push `(sqrt(*xmax*) = sqrt(,*xmax*) = ,y) *results*)
	     (push '(test of error returns) *results*)
	     (setq x 0.0)
	     (push `(sqrt will be called with an argument of ,x this
			   should not trigger an error) *results*)
	     (setq y (sqrt x))
	     (push `(sqrt returned the value ,y) *results*)
 	     (setq x -1.0)
 	     (push `(sqrt will be called with an argument of ,x this
 			   should trigger an error) *results*)
 	     (setq y (car (errset (sqrt x))))
 	     (push `(sqrt returned the value ,y) *results*)
	     (push '(this concludes the tests) *results*)
	     t))
;;; atan Test Routine

(defmacro atan2 (x y) `(atan ,x ,y))

(defun atan-test ()
       (declare (flonum beta albeta ait a b betap del em expon ob32 ran r6 r7
			xn w x xl xsq x1 y z zz)
		(fixnum n i k1 k2 k3 ii i1 j))
       (machar)
       (let* ((beta (float *ibeta*))
	      (albeta (log beta))(ait (float *it*))
	      (a -0.0625)(b (minus a))(ob32 (*$ b 0.5))
	      (n 8000.)(xn (float n))(x1 0.0)(xl 0.0)(expon 0.0)
 	      (r7 0.0)(r6 0.0)(x 0.0)(z 0.0)(xsq 0.0)(em 0.0)
	      (k1 0)(k3 0)(del 0.0)(sum 0.0)(zz 0.0)
	      (i1 0)(y 0.0)(w 0.0)(k2 0)(betap 0.0))
	     (do ((j 1 (1+ j)))
		 ((> j 4))
		 (setq k1 0)(setq x1 0.0)(setq r7 0.0)
		 (setq k3 0)(setq r6 0.0)(setq del (//$ (-$ b a) xn))
		 (setq xl a)
		 (do ((i 1 (1+ i)))
		     ((> i n))
		     (setq x (+$ (*$ del (ran i1))
				 xl))
		     (cond ((= j 2)
			    (setq x (*$ (-$ (+$ 1.0 (*$ x a)) 1.0) 16.0))))
		     (setq z (atan x 1.0))
		     (cond ((= j 1)
			    (setq xsq (*$ x x))
			    (setq em 17.0)
			    (setq sum (//$ xsq em))
			    (do ((ii 1 (1+ ii)))
				((> ii 7))
				(setq em (-$ em 2.0))
				(setq sum (*$ (-$ (//$ 1.0 em) sum) xsq)))
			    (setq sum (*$ (minus x) sum))
			    (setq zz (+$ x sum))
			    (cond ((zerop *irnd*) (setq zz (+$ zz (+$ sum sum))))))
			   (t (cond ((= j 2)
				     (setq y (-$ x .0625))
				     (setq y (//$ y (+$ 1.0 (*$ x a))))
				     (setq zz (+$ (-$ (atan y 1.0)
						      8.11900402651526021e-5)
						  ob32))
				     (setq zz (+$ zz ob32)))
				    (t (setq z (+$ z z))
				       (setq
					y (//$ x (*$ (+$ 0.5 (*$ x 0.5))
						     (+$ (-$ 0.5 x) 0.5))))
				       (setq zz (atan y 1.0))))))
		     (setq w 1.0)
		     (cond ((not (zerop z))
			    (setq w (//$ (-$ z zz) z))))
		     (cond ((> w 0.0)
			    (setq k1 (1+ k1))))
		     (cond ((< w 0.0)
			    (setq k3 (1+ k3))))
		     (setq w (abs w))
		     (cond ((< r6 w)
			    (setq r6 w)
			    (setq x1 x)))
		     (setq r7 (+$ r7
				  (*$ w w)))
		     (setq xl (+$ xl del)))
	     (setq k2 (- (- n k3) k1))
	     (setq r7 (sqrt (//$ r7 xn)))
	     (cond ((= j 1)
		    (push `(Test of atan ( x ) vs truncated
				 taylor series)
			  *results*)))
	     (cond ((= j 2)
		    (push `(Test of atan  ( x )
				 vs atan ( 1 // 16.) +
				 atan((x - 1 // 16.)//(1 + x // 16.)))
			  *results*)))
	     (cond ((= j 3)
		    (push `(test of 2 * atan  ( x ) vs atan
				 (2x // (1 - x * x)))
			  *results*)))
	     (push `(,n random arguments were tested from the
			interval ( ,a ,b))
		   *results*)
	     (push `(atan ( x )  was larger ,k1 times)
		   *results*)
	     (push `(it agreed ,k2 times) *results*)
	     (push `(it was smaller ,k3 times) *results*)
	     (push `(there are ,*it* significant base ,*ibeta* digits in
			   a floating-point number) *results*)
	     (setq w -999.0)
	     (cond ((not (zerop r6))
		    (setq w (//$ (log (abs r6)) albeta))))
	     (push `(the maximum relative error of ,r6 = ,*ibeta*
			 ↑ ,w occurred for x = ,x1)
		   *results*)
	     (setq w (max (+$ ait w) 0.0))
	     (push `(the estimated loss of base ,*ibeta* significant digits
			 is ,w) *results*)
	     (setq w -999.0)
	     (cond ((not (zerop r7))
		    (setq w (//$ (log (abs r7)) albeta))))
	     (push `(the root mean square relative error was ,r7
			 = ,*ibeta* ↑ ,w) *results*)
	     (setq w (max (+$ ait w) 0.0))
	     (push `(the estimated loss of base ,*ibeta* significant digits
			 is ,w) *results*)
	     (setq a b)
	     (cond ((= j 1) (setq b (-$ 2.0 (sqrt 3.0)))))
	     (cond ((= j 2) (setq b (-$ (sqrt 2.0) 1.0))))
	     (cond ((= j 3) (setq b 1.0))))
       (push `(special tests) *results*)
       (push `(the identity: atan (-x) = -atan (x) will be tested)
	     *results*)
       (push `(x : f (x) + f (-x)) *results*)
       (setq a 5.0)
       (do ((i 1 (1+ i)))
	   ((> i 5))
	   (setq x (*$ (ran i1) a))
	   (setq z (+$ (atan x 1.0)
		       (atan (minus x) 1.0)))
	   (push `(,x : ,z) *results*))
       (push `(the identity atan (x) = x for x small will be tested)
	     *results*)
       (push `(x : x - f ( x)) *results*)
       (setq betap (expt beta *it*))
       (setq x (//$ (ran i1) betap))
       (do ((i 1 (1+ i)))
	   ((> i 5))
	   (setq z (-$ x (atan x 1.0)))
	   (push `(,x : ,z) *results*)
	   (setq x (//$ x beta)))
       (push `(the identity atan (x // y) = atan2 (x y) will
		   be tested) *results*)
       (push `(the first column of results should be 0 and the
		   second should be +-π) *results*)
       (push `(x : y : f1 ( x // y) - f2 (x y) : f1 (x // y) - f2 (x // -y))
	     *results*)
       (setq a -2.0)
       (setq  b 4.0)
       (do ((i 1 (1+ i)))
	   ((> i 5))
	   (setq x (+$ (*$ (ran i1) b) a))
	   (setq y (ran i1))
	   (setq w (minus y))
	   (setq z (-$ (atan (//$ x y) 1.0) (atan2 x y)))
	   (setq zz (-$ (atan (//$ x w) 1.0) (atan2 x w)))
	   (push `(,x : ,y : ,z : ,zz) *results*))
       (push `(test of very small argument) *results*)
       (setq expon (*$ (float *minexp*) 0.75))
       (setq x (expt beta expon))
       (setq y (atan x 1.0))
       (push `(atan ( ,x ) = ,y) *results*)
       (push `(test of error returns) *results*)
       (push `(atan will be called with the argument ,*xmax*) *results*)
       (push `(this should not trigger an error message) *results*)
       (setq z (car (errset (atan *xmax* 1.0))))
       (push `(atan ( ,*xmax* ) = ,z) *results*)
       (setq x 1.0 y 0.0)
       (push `(atan2 will be called with the arguments ,x ,y) *results*)
       (push `(this should not trigger an error message) *results*)
       (setq z (car (errset (atan2 x y))))
       (push `(atan2 ( ,x ,y) = ,z) *results*)
       (push `(atan2 will be called with the arguments ,*xmin*
		       ,*xmax*) *results*)
       (push `(this should not trigger an error message) *results*)
       (setq z (car (errset (atan2 *xmin* *xmax*))))
       (push `(atan2 ( ,*xmin* ,*xmax*) = ,z) *results*)
       (push `(atan2 will be called with the arguments ,*xmax*
		       ,*xmin*) *results*)
       (push `(this should not trigger an error message) *results*)
       (setq z (car (errset (atan2 *xmax* *xmin*))))
       (push `(atan2 ( ,*xmax* ,*xmin*) = ,z) *results*)
       (setq x 0.0)
       (push `(atan2 will be called with the arguments ,x ,y)
 	     *results*)
       (push `(this should trigger an error message) *results*)
       (setq z (car (errset (atan2 x y))))
       (push `(atan2 ( ,x ,y) = ,z) *results*)
       (push `(This concludes the tests) *results*)
       t))
;;; Testing

(defun show-results ()
       (mapc #'print (reverse *results*)) t)

(include "timer.lsp[tim,lsp]")

(timer timit1 (progn (sqrt-test)(atan-test)))

(defun timit ()
       (timit1)
       (show-results))